clear all ; close all ;


ggg = dir('*_Pleth_PTT.txt') ;
addpath('/Users/hautiengwu/Desktop/dsSTFT_code/CodeGeneric') ;

Hz = 4 ;
basicTF.win = 1000 ; %100*4*10 ; %4096;
basicTF.hop = 10; %441;
basicTF.fs = Hz;
basicTF.fr = 2e-4;
basicTF.feat = 'SST11';
advTF.num_tap = 1;%num_tap(cc);
advTF.win_type = 'Gauss'; %{'Gaussian','Thomson','multipeak','SWCE'}; %}
advTF.Smo = 1;
advTF.Rej = 0;
advTF.ths = 1E-9;
advTF.HighFreq = 1/50;
advTF.LowFreq = 0.005/50;
advTF.lpc = 0;
cepR.g = 0.2 ;
cepR.Tc=0;
P.num_s = 1;
P.num_c = 1;





fprintf(['Case NO = ',num2str(length(ggg)),'\n']) ;
LEN = zeros(length(ggg),1) ;


for caseNO = 1:length(ggg)        
    [x0] = dlmread(ggg(caseNO).name) ;
    LEN(caseNO) = size(x0,1) ;
    
    sig = x0(:,2)-median(x0(:,2)) ;
    trend = zeros(size(sig)) ;
    for jj =1 :length(trend)
        idx = max(1, jj-100*4) : min(length(sig), jj+100*4) ;
        trend(jj) = median(sig(idx)) ;
    end
    Jidx = find(abs(sig-trend)<25) ;
    sig2 = interp1(Jidx, sig(Jidx)-trend(Jidx), [1:length(sig)], 'linear', 'extrap') ;

    ALLsig2{caseNO} = sig2 ;
    ALLsig{caseNO} = sig-trend ;
    ALLtrend{caseNO} = trend ;
end



%% plot signal
figure
for ii = 1:length(ggg) ; 
    plot([1:4:length(ALLsig{ii})]/4, ALLsig{ii}(1:4:end)+ALLtrend{ii}(1:4:end)+(ii-1)*30) ; hold on; 
end
axis([-inf inf -30 700]) ; set(gca,'fontsize', 24) ;
xlabel('Time (sec)') ; ylabel('a.u.') ;
export_fig(['AllSignalFigure'],'-transparent','-pdf','-nocrop');



%% plot spectrum
figure
for ii=1:length(ggg)
    xi = 2*[1:(LEN(ii)-1)/2]./((LEN(ii)-1)/2);
    xhat = fft(ALLsig2{ii}) ; 
    [a,b]=max(abs(xhat(2:(LEN(ii)-1)/2)));
    PEAK(ii) = 1./xi(b(1)) ; 
    plot(xi, 1e-4*(3*abs(xhat(2:(LEN(ii)-1)/2+1))+(ii-1)*.5e4)) ; hold on; 
end 
plot([0.01 0.01], [0 12], 'k')
axis([0 0.12 -inf inf] );
set(gca,'fontsize', 24) ;
xlabel('Frequency (Hz)'); ylabel('a.u.') ;
export_fig(['AllSignalFigurePS'],'-transparent','-pdf','-nocrop');




%% dsSST

MEANdeShapeSTFT = zeros(400, 62) ;
MEANdeShapeSST = zeros(400, 62) ;

for caseNO = 1:length(ggg)
    
    disp(caseNO) ;
    
    [x0] = dlmread(ggg(caseNO).name) ;
    LEN(caseNO) = size(x0,1) ;
    
    sig = x0(:,2)-median(x0(:,2)) ;
    trend = zeros(size(sig)) ;
    for jj =1 :length(trend)
        idx = max(1, jj-100*4) : min(length(sig), jj+100*4) ;
        trend(jj) = median(sig(idx)) ;
    end
    sig0 = sig-trend ;
    sig = sig-trend ;
    
    for jj = 1:length(sig)
        idx = max(1, jj-60*4) : min(length(sig), jj+60*4) ;
        if sig(jj)>quantile(sig(idx), .92) * 2 ;
            sig(jj) = median(sig) ;
        elseif sig(jj)<quantile(sig(idx), .05)- 10;
            sig(jj) = median(sig) ;
        end
    end
    
    x = sig ; x(find(x>50)) = 50 ;
    t = [1:length(sig)]/ Hz ;
    disp(length(x)) ;
    
    if length(x)<1201
        fprintf('TOO Short signal\n') ;
    else
        
        time1 = tic ;
        [~, ~, ~, deShapeSTFT, deShapeSST, tfrsq, tfrtic] = CFPH(x+1, basicTF, advTF, cepR, P);
        toc(time1) 
    
        time_stamp = basicTF.hop/basicTF.fs;
    
        figure ; 
        imageSQ(0:time_stamp:t(end), tfrtic*basicTF.fs, abs(deShapeSST), 0.995); axis xy; colormap(1-gray);
        xlabel('time (sec)'); ylabel('Frequency (Hz)'); set(gca,'fontsize',28) ; 
        set(gca,'ytick',[0.01:0.01:0.08]) ; %title('MEAN SST') ;
        yyaxis right
        plot(t, x0(:,2),'b'); axis([-inf inf 100 300]);
        ylabel('PTT (ms)') ;
        export_fig(['PTT_deShapeSTFT_case',num2str(caseNO)],'-transparent','-pdf','-nocrop');
        close 
    
        MEANdeShapeSTFT = MEANdeShapeSTFT + abs(deShapeSTFT) ;    
        MEANdeShapeSST = MEANdeShapeSST + abs(deShapeSST) ;
    end
    
    clear deShapeSTFT; clear deShapeSST ; clear x ;
end




% subplot(221) ; hold off;
% imageSQ(0:time_stamp:t(end), tfrtic*basicTF.fs, MAXdeShapeSTFT, 0.995); axis xy; colormap(1-gray);
% xlabel('time (sec)'); ylabel('Frequency (Hz)'); set(gca,'fontsize',28) ; title('Max STFT') ;
% 
% subplot(222) ; hold off;
% imageSQ(0:time_stamp:t(end), tfrtic*basicTF.fs, MAXdeShapeSST, 0.995); axis xy; colormap(1-gray);
% xlabel('time (sec)'); ylabel('Frequency (Hz)'); set(gca,'fontsize',28) ; title('Max SST') ;
% 
% subplot(223) ; hold off;
% imageSQ(0:time_stamp:t(end), tfrtic*basicTF.fs, MEANdeShapeSTFT, 0.995); axis xy; colormap(1-gray);
% xlabel('time (sec)'); ylabel('Frequency (Hz)'); set(gca,'fontsize',28) ; title('Mean STFT') ;

%subplot(224) ; hold off;
figure ; imageSQ(0:time_stamp:t(end), tfrtic*basicTF.fs, MEANdeShapeSST, 0.995); axis xy; colormap(1-gray);
xlabel('time (sec)'); ylabel('Frequency (Hz)'); set(gca,'fontsize',28) ; 
set(gca,'ytick',[0.01:0.01:0.08])
%title('MEAN SST') ;

%yyaxis right
%plot(t, x0(:,2),'b'); axis([-inf inf 100 300]);

export_fig('PTT_deShapeSTFT','-transparent','-pdf','-nocrop');

